SUBROUTINE qcim_dalton_initial

   USE qcim_main
   IMPLICIT NONE
   INTEGER :: nbas, natoms, nshells
   INTEGER, PARAMETER :: IO2 = 777

   OPEN(UNIT=IO2, FILE='basis.qcim', STATUS='OLD',  &
        ACCESS='SEQUENTIAL', FORM='FORMATTED')
   READ(IO2,'(I8)') natoms
   READ(IO2,'(I8)') nshells
   READ(IO2,'(I8)') nbas
   CLOSE(IO2)
   CALL qcim_initial(nbas)

END SUBROUTINE qcim_dalton_initial

!-------------------------------------------------------------------------------

SUBROUTINE qcim_dalton_final
   USE qcim_main
   IMPLICIT NONE
   CALL qcim_final
END SUBROUTINE qcim_dalton_final

!-------------------------------------------------------------------------------

SUBROUTINE qcim_dalton_run(n,nbas,DMAT,QMAT)

   USE qcim_global_consts
   USE qcim_main, ONLY: qcim_get_Q_matrix
   USE fmm_main, ONLY: fmm_get_J_matrix

   IMPLICIT NONE
   INTEGER, INTENT(IN)  :: n, nbas
   REAL(8), INTENT(IN)  :: DMAT(n)
   REAL(8), INTENT(OUT) :: QMAT(n)

   REAL(8) :: dens(nbas,nbas), Q_matrix(nbas,nbas)
   INTEGER :: p,q,k

   ! Factor of 2 for off-diagonal density elements assumed to be
   ! included in folded input format of DMAT. 
   ! We also keep this factor in the unpacked square form since we only
   ! loop over upper triangular pairs in general.
   Q_matrix(:,:) = zero
   dens(:,:) = zero
   k = 0
   DO q = 1, nbas
      DO p = 1, q
         k = k + 1
         IF (k>n) STOP 'k>n in qcim_dalton_run'
         dens(p,q) = DMAT(k)
         dens(q,p) = DMAT(k)
      END DO
   END DO

   print *, "Getting QCIM update via Laplace equation..."
   CALL qcim_get_Q_matrix(nbas,dens,Q_matrix)

   ! print *, "Getting QCIM update via FMM engine..."
   ! (Don't forget to initialize fmm with 'MD4_FMM')
   ! CALL fmm_get_J_matrix('ONE_EL',dens,Q_matrix)

   ! Dalton expects a negative correction (nuclear-attraction term) due
   ! to external positive charges so we need a minus sign on Q_matrix
   k = 0
   DO q = 1, nbas
      DO p = 1, q
         k = k + 1
         IF (k>n) STOP 'k>n in qcim_dalton_run'
         QMAT(k) = -Q_matrix(p,q)
      END DO
   END DO


END SUBROUTINE qcim_dalton_run

!-------------------------------------------------------------------------------
